<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<title>Voro++: unitcell.cc Source File</title>

<link href="tabs.css" rel="stylesheet" type="text/css"/>
<link href="doxygen.css" rel="stylesheet" type="text/css" />



</head>
<body>
<div id="top"><!-- do not remove this div! -->


<div id="titlearea">
<table cellspacing="0" cellpadding="0">
 <tbody>
 <tr style="height: 56px;">
  
  
  <td style="padding-left: 0.5em;">
   <div id="projectname">Voro++
   
   </div>
   
  </td>
  
  
  
 </tr>
 </tbody>
</table>
</div>

<!-- Generated by Doxygen 1.7.5.1 -->
<script type="text/javascript" src="dynsections.js"></script>
  <div id="navrow1" class="tabs">
    <ul class="tablist">
      <li><a href="index.html"><span>Main&#160;Page</span></a></li>
      <li><a href="annotated.html"><span>Data&#160;Structures</span></a></li>
      <li class="current"><a href="files.html"><span>Files</span></a></li>
    </ul>
  </div>
  <div id="navrow2" class="tabs2">
    <ul class="tablist">
      <li><a href="files.html"><span>File&#160;List</span></a></li>
      <li><a href="globals.html"><span>Globals</span></a></li>
    </ul>
  </div>
<div class="header">
  <div class="headertitle">
<div class="title">unitcell.cc</div>  </div>
</div>
<div class="contents">
<a href="unitcell_8cc.html">Go to the documentation of this file.</a><div class="fragment"><pre class="fragment"><a name="l00001"></a>00001 <span class="comment">// Voro++, a 3D cell-based Voronoi library</span>
<a name="l00002"></a>00002 <span class="comment">//</span>
<a name="l00003"></a>00003 <span class="comment">// Author   : Chris H. Rycroft (LBL / UC Berkeley)</span>
<a name="l00004"></a>00004 <span class="comment">// Email    : chr@alum.mit.edu</span>
<a name="l00005"></a>00005 <span class="comment">// Date     : August 30th 2011</span>
<a name="l00006"></a>00006 <span class="comment"></span>
<a name="l00007"></a>00007 <span class="comment">/** \file unitcell.cc</span>
<a name="l00008"></a>00008 <span class="comment"> * \brief Function implementations for the unitcell class. */</span>
<a name="l00009"></a>00009 
<a name="l00010"></a>00010 <span class="preprocessor">#include &lt;queue&gt;</span>
<a name="l00011"></a>00011 <span class="keyword">using namespace </span>std;
<a name="l00012"></a>00012 
<a name="l00013"></a>00013 <span class="preprocessor">#include &quot;<a class="code" href="unitcell_8hh.html" title="Header file for the unitcell class.">unitcell.hh</a>&quot;</span>
<a name="l00014"></a>00014 <span class="preprocessor">#include &quot;<a class="code" href="cell_8hh.html" title="Header file for the voronoicell and related classes.">cell.hh</a>&quot;</span>
<a name="l00015"></a>00015 
<a name="l00016"></a>00016 <span class="keyword">namespace </span>voro {
<a name="l00017"></a>00017 <span class="comment"></span>
<a name="l00018"></a>00018 <span class="comment">/** Initializes the unit cell class for a particular non-orthogonal periodic</span>
<a name="l00019"></a>00019 <span class="comment"> * geometry, corresponding to a parallelepiped with sides given by three</span>
<a name="l00020"></a>00020 <span class="comment"> * vectors. The class constructs the unit Voronoi cell corresponding to this</span>
<a name="l00021"></a>00021 <span class="comment"> * geometry.</span>
<a name="l00022"></a>00022 <span class="comment"> * \param[in] (bx_) The x coordinate of the first unit vector.</span>
<a name="l00023"></a>00023 <span class="comment"> * \param[in] (bxy_,by_) The x and y coordinates of the second unit vector.</span>
<a name="l00024"></a>00024 <span class="comment"> * \param[in] (bxz_,byz_,bz_) The x, y, and z coordinates of the third unit</span>
<a name="l00025"></a>00025 <span class="comment"> *                            vector. */</span>
<a name="l00026"></a><a class="code" href="classvoro_1_1unitcell.html#abbeb66edf9e36476cf47bee76ce8d9fc">00026</a> unitcell::unitcell(<span class="keywordtype">double</span> bx_,<span class="keywordtype">double</span> bxy_,<span class="keywordtype">double</span> by_,<span class="keywordtype">double</span> bxz_,<span class="keywordtype">double</span> byz_,<span class="keywordtype">double</span> bz_)
<a name="l00027"></a>00027         : bx(bx_), bxy(bxy_), by(by_), bxz(bxz_), byz(byz_), bz(bz_) {
<a name="l00028"></a>00028         <span class="keywordtype">int</span> i,j,l=1;
<a name="l00029"></a>00029 
<a name="l00030"></a>00030         <span class="comment">// Initialize the Voronoi cell to be a very large rectangular box</span>
<a name="l00031"></a>00031         <span class="keyword">const</span> <span class="keywordtype">double</span> ucx=max_unit_voro_shells*<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>,ucy=max_unit_voro_shells*<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>,ucz=max_unit_voro_shells*<a class="code" href="classvoro_1_1unitcell.html#a30bf0a4a6bff808301aa424683eb0481">bz</a>;
<a name="l00032"></a>00032         <a class="code" href="classvoro_1_1unitcell.html#a0844f70778cc7b19024ed28a2c7ad02c">unit_voro</a>.<a class="code" href="classvoro_1_1voronoicell.html#a573257366d690af4fd8646b4c431ca84">init</a>(-ucx,ucx,-ucy,ucy,-ucz,ucz);
<a name="l00033"></a>00033 
<a name="l00034"></a>00034         <span class="comment">// Repeatedly cut the cell by shells of periodic image particles</span>
<a name="l00035"></a>00035         <span class="keywordflow">while</span>(l&lt;2*max_unit_voro_shells) {
<a name="l00036"></a>00036 
<a name="l00037"></a>00037                 <span class="comment">// Check to see if any of the planes from the current shell</span>
<a name="l00038"></a>00038                 <span class="comment">// will cut the cell</span>
<a name="l00039"></a>00039                 <span class="keywordflow">if</span>(unit_voro_intersect(l)) {
<a name="l00040"></a>00040 
<a name="l00041"></a>00041                         <span class="comment">// If they do, apply the plane cuts from the current</span>
<a name="l00042"></a>00042                         <span class="comment">// shell</span>
<a name="l00043"></a>00043                         unit_voro_apply(l,0,0);
<a name="l00044"></a>00044                         <span class="keywordflow">for</span>(i=1;i&lt;l;i++) {
<a name="l00045"></a>00045                                 unit_voro_apply(l,i,0);
<a name="l00046"></a>00046                                 unit_voro_apply(-l,i,0);
<a name="l00047"></a>00047                         }
<a name="l00048"></a>00048                         <span class="keywordflow">for</span>(i=-l;i&lt;=l;i++) unit_voro_apply(i,l,0);
<a name="l00049"></a>00049                         <span class="keywordflow">for</span>(i=1;i&lt;l;i++) <span class="keywordflow">for</span>(j=-l+1;j&lt;=l;j++) {
<a name="l00050"></a>00050                                 unit_voro_apply(l,j,i);
<a name="l00051"></a>00051                                 unit_voro_apply(-j,l,i);
<a name="l00052"></a>00052                                 unit_voro_apply(-l,-j,i);
<a name="l00053"></a>00053                                 unit_voro_apply(j,-l,i);
<a name="l00054"></a>00054                         }
<a name="l00055"></a>00055                         <span class="keywordflow">for</span>(i=-l;i&lt;=l;i++) <span class="keywordflow">for</span>(j=-l;j&lt;=l;j++) unit_voro_apply(i,j,l);
<a name="l00056"></a>00056                 } <span class="keywordflow">else</span> {
<a name="l00057"></a>00057 
<a name="l00058"></a>00058                         <span class="comment">// Calculate a bound on the maximum y and z coordinates</span>
<a name="l00059"></a>00059                         <span class="comment">// that could possibly cut the cell. This is based upon</span>
<a name="l00060"></a>00060                         <span class="comment">// a geometric result that particles with z&gt;l can&#39;t cut</span>
<a name="l00061"></a>00061                         <span class="comment">// a cell lying within the paraboloid</span>
<a name="l00062"></a>00062                         <span class="comment">// z&lt;=(l*l-x*x-y*y)/(2*l). It is always a tighter bound</span>
<a name="l00063"></a>00063                         <span class="comment">// than the one based on computing the maximum radius</span>
<a name="l00064"></a>00064                         <span class="comment">// of a Voronoi cell vertex.</span>
<a name="l00065"></a>00065                         <a class="code" href="classvoro_1_1unitcell.html#a3f7cebe966d72935971c508fc64902f0">max_uv_y</a>=<a class="code" href="classvoro_1_1unitcell.html#a3ac0618230ffa0a382469530f97844e6">max_uv_z</a>=0;
<a name="l00066"></a>00066                         <span class="keywordtype">double</span> y,z,q,*pts=<a class="code" href="classvoro_1_1unitcell.html#a0844f70778cc7b19024ed28a2c7ad02c">unit_voro</a>.<a class="code" href="classvoro_1_1voronoicell__base.html#aa4e3f22d4ff450097571e95caeda1b3d">pts</a>,*pp=pts;
<a name="l00067"></a>00067                         <span class="keywordflow">while</span>(pp&lt;pts+3*<a class="code" href="classvoro_1_1unitcell.html#a0844f70778cc7b19024ed28a2c7ad02c">unit_voro</a>.<a class="code" href="classvoro_1_1voronoicell__base.html#a7e69469f95464afbeb9feec927507243">p</a>) {
<a name="l00068"></a>00068                                 q=*(pp++);y=*(pp++);z=*(pp++);q=sqrt(q*q+y*y+z*z);
<a name="l00069"></a>00069                                 <span class="keywordflow">if</span>(y+q&gt;<a class="code" href="classvoro_1_1unitcell.html#a3f7cebe966d72935971c508fc64902f0">max_uv_y</a>) <a class="code" href="classvoro_1_1unitcell.html#a3f7cebe966d72935971c508fc64902f0">max_uv_y</a>=y+q;
<a name="l00070"></a>00070                                 <span class="keywordflow">if</span>(z+q&gt;<a class="code" href="classvoro_1_1unitcell.html#a3ac0618230ffa0a382469530f97844e6">max_uv_z</a>) <a class="code" href="classvoro_1_1unitcell.html#a3ac0618230ffa0a382469530f97844e6">max_uv_z</a>=z+q;
<a name="l00071"></a>00071                         }
<a name="l00072"></a>00072                         <a class="code" href="classvoro_1_1unitcell.html#a3ac0618230ffa0a382469530f97844e6">max_uv_z</a>*=0.5;
<a name="l00073"></a>00073                         <a class="code" href="classvoro_1_1unitcell.html#a3f7cebe966d72935971c508fc64902f0">max_uv_y</a>*=0.5;
<a name="l00074"></a>00074                         <span class="keywordflow">return</span>;
<a name="l00075"></a>00075                 }
<a name="l00076"></a>00076                 l++;
<a name="l00077"></a>00077         }
<a name="l00078"></a>00078 
<a name="l00079"></a>00079         <span class="comment">// If the routine makes it here, then the unit cell still hasn&#39;t been</span>
<a name="l00080"></a>00080         <span class="comment">// completely bounded by the plane cuts. Give the memory error code,</span>
<a name="l00081"></a>00081         <span class="comment">// because this is mainly a case of hitting a safe limit, than any</span>
<a name="l00082"></a>00082         <span class="comment">// inherent problem.</span>
<a name="l00083"></a>00083         voro_fatal_error(<span class="stringliteral">&quot;Periodic cell computation failed&quot;</span>,<a class="code" href="config_8hh.html#a5dc0616f8a67ae3d1c2ba8a3dcf5612b">VOROPP_MEMORY_ERROR</a>);
<a name="l00084"></a>00084 }
<a name="l00085"></a>00085 <span class="comment"></span>
<a name="l00086"></a>00086 <span class="comment">/** Applies a pair of opposing plane cuts from a periodic image point</span>
<a name="l00087"></a>00087 <span class="comment"> * to the unit Voronoi cell.</span>
<a name="l00088"></a>00088 <span class="comment"> * \param[in] (i,j,k) the index of the periodic image to consider. */</span>
<a name="l00089"></a>00089 <span class="keyword">inline</span> <span class="keywordtype">void</span> unitcell::unit_voro_apply(<span class="keywordtype">int</span> i,<span class="keywordtype">int</span> j,<span class="keywordtype">int</span> k) {
<a name="l00090"></a>00090         <span class="keywordtype">double</span> x=i*<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+j*<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>+k*<a class="code" href="classvoro_1_1unitcell.html#a000adf05e666f0fb4f57dbe8ca5b51db">bxz</a>,y=j*<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>+k*<a class="code" href="classvoro_1_1unitcell.html#a434c1a7cd8ba84d4ad09845a1a0ddb5b">byz</a>,z=k*<a class="code" href="classvoro_1_1unitcell.html#a30bf0a4a6bff808301aa424683eb0481">bz</a>;
<a name="l00091"></a>00091         <a class="code" href="classvoro_1_1unitcell.html#a0844f70778cc7b19024ed28a2c7ad02c">unit_voro</a>.<a class="code" href="classvoro_1_1voronoicell.html#a6e95394d82fd835c75a77462c872594e">plane</a>(x,y,z);
<a name="l00092"></a>00092         <a class="code" href="classvoro_1_1unitcell.html#a0844f70778cc7b19024ed28a2c7ad02c">unit_voro</a>.<a class="code" href="classvoro_1_1voronoicell.html#a6e95394d82fd835c75a77462c872594e">plane</a>(-x,-y,-z);
<a name="l00093"></a>00093 }
<a name="l00094"></a>00094 <span class="comment"></span>
<a name="l00095"></a>00095 <span class="comment">/** Calculates whether the unit Voronoi cell intersects a given periodic image</span>
<a name="l00096"></a>00096 <span class="comment"> * of the domain.</span>
<a name="l00097"></a>00097 <span class="comment"> * \param[in] (dx,dy,dz) the displacement of the periodic image.</span>
<a name="l00098"></a>00098 <span class="comment"> * \param[out] vol the proportion of the unit cell volume within this image,</span>
<a name="l00099"></a>00099 <span class="comment"> *                 only computed in the case that the two intersect.</span>
<a name="l00100"></a>00100 <span class="comment"> * \return True if they intersect, false otherwise. */</span>
<a name="l00101"></a><a class="code" href="classvoro_1_1unitcell.html#acc445d2a964296dfca277e03fdc63ebd">00101</a> <span class="keywordtype">bool</span> <a class="code" href="classvoro_1_1unitcell.html#acc445d2a964296dfca277e03fdc63ebd">unitcell::intersects_image</a>(<span class="keywordtype">double</span> dx,<span class="keywordtype">double</span> dy,<span class="keywordtype">double</span> dz,<span class="keywordtype">double</span> &amp;vol) {
<a name="l00102"></a>00102         <span class="keyword">const</span> <span class="keywordtype">double</span> bxinv=1/<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>,byinv=1/<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>,bzinv=1/<a class="code" href="classvoro_1_1unitcell.html#a30bf0a4a6bff808301aa424683eb0481">bz</a>,ivol=bxinv*byinv*bzinv;
<a name="l00103"></a>00103         <a class="code" href="classvoro_1_1voronoicell.html" title="Extension of the voronoicell_base class to represent a Voronoi cell without neighbor information...">voronoicell</a> c;
<a name="l00104"></a>00104         c=<a class="code" href="classvoro_1_1unitcell.html#a0844f70778cc7b19024ed28a2c7ad02c">unit_voro</a>;
<a name="l00105"></a>00105         dx*=2;dy*=2;dz*=2;
<a name="l00106"></a>00106         <span class="keywordflow">if</span>(!c.<a class="code" href="classvoro_1_1voronoicell.html#a6e95394d82fd835c75a77462c872594e">plane</a>(0,0,bzinv,dz+1)) <span class="keywordflow">return</span> <span class="keyword">false</span>;
<a name="l00107"></a>00107         <span class="keywordflow">if</span>(!c.<a class="code" href="classvoro_1_1voronoicell.html#a6e95394d82fd835c75a77462c872594e">plane</a>(0,0,-bzinv,-dz+1)) <span class="keywordflow">return</span> <span class="keyword">false</span>;
<a name="l00108"></a>00108         <span class="keywordflow">if</span>(!c.<a class="code" href="classvoro_1_1voronoicell.html#a6e95394d82fd835c75a77462c872594e">plane</a>(0,byinv,-byz*byinv*bzinv,dy+1)) <span class="keywordflow">return</span> <span class="keyword">false</span>;
<a name="l00109"></a>00109         <span class="keywordflow">if</span>(!c.<a class="code" href="classvoro_1_1voronoicell.html#a6e95394d82fd835c75a77462c872594e">plane</a>(0,-byinv,byz*byinv*bzinv,-dy+1)) <span class="keywordflow">return</span> <span class="keyword">false</span>;
<a name="l00110"></a>00110         <span class="keywordflow">if</span>(!c.<a class="code" href="classvoro_1_1voronoicell.html#a6e95394d82fd835c75a77462c872594e">plane</a>(bxinv,-<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>*bxinv*byinv,(<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>*byz-<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>*bxz)*ivol,dx+1)) <span class="keywordflow">return</span> <span class="keyword">false</span>;
<a name="l00111"></a>00111         <span class="keywordflow">if</span>(!c.<a class="code" href="classvoro_1_1voronoicell.html#a6e95394d82fd835c75a77462c872594e">plane</a>(-bxinv,<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>*bxinv*byinv,(-<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>*byz+<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>*bxz)*ivol,-dx+1)) <span class="keywordflow">return</span> <span class="keyword">false</span>;
<a name="l00112"></a>00112         vol=c.<a class="code" href="classvoro_1_1voronoicell__base.html#a92f69986a6f78e7b8b6d9f5778ed99ba">volume</a>()*ivol;
<a name="l00113"></a>00113         <span class="keywordflow">return</span> <span class="keyword">true</span>;
<a name="l00114"></a>00114 }
<a name="l00115"></a>00115 <span class="comment"></span>
<a name="l00116"></a>00116 <span class="comment">/** Computes a list of periodic domain images that intersect the unit Voronoi cell.</span>
<a name="l00117"></a>00117 <span class="comment"> * \param[out] vi a vector containing triplets (i,j,k) corresponding to domain</span>
<a name="l00118"></a>00118 <span class="comment"> *                images that intersect the unit Voronoi cell, when it is</span>
<a name="l00119"></a>00119 <span class="comment"> *                centered in the middle of the primary domain.</span>
<a name="l00120"></a>00120 <span class="comment"> * \param[out] vd a vector containing the fraction of the Voronoi cell volume</span>
<a name="l00121"></a>00121 <span class="comment"> *                within each corresponding image listed in vi. */</span>
<a name="l00122"></a><a class="code" href="classvoro_1_1unitcell.html#abeb3c87fffa74c969bd180a5f1f0ccb2">00122</a> <span class="keywordtype">void</span> <a class="code" href="classvoro_1_1unitcell.html#abeb3c87fffa74c969bd180a5f1f0ccb2">unitcell::images</a>(vector&lt;int&gt; &amp;vi,vector&lt;double&gt; &amp;vd) {
<a name="l00123"></a>00123         <span class="keyword">const</span> <span class="keywordtype">int</span> ms2=max_unit_voro_shells*2+1,mss=ms2*ms2*ms2;
<a name="l00124"></a>00124         <span class="keywordtype">bool</span> *a=<span class="keyword">new</span> <span class="keywordtype">bool</span>[mss],*ac=a+max_unit_voro_shells*(1+ms2*(1+ms2)),*ap=a;
<a name="l00125"></a>00125         <span class="keywordtype">int</span> i,j,k;
<a name="l00126"></a>00126         <span class="keywordtype">double</span> vol;
<a name="l00127"></a>00127 
<a name="l00128"></a>00128         <span class="comment">// Initialize mask</span>
<a name="l00129"></a>00129         <span class="keywordflow">while</span>(ap&lt;ac) *(ap++)=<span class="keyword">true</span>;
<a name="l00130"></a>00130         *(ap++)=<span class="keyword">false</span>;
<a name="l00131"></a>00131         <span class="keywordflow">while</span>(ap&lt;a+mss) *(ap++)=<span class="keyword">true</span>;
<a name="l00132"></a>00132 
<a name="l00133"></a>00133         <span class="comment">// Set up the queue and add (0,0,0) image to it</span>
<a name="l00134"></a>00134         queue&lt;int&gt; q;
<a name="l00135"></a>00135         q.push(0);q.push(0);q.push(0);
<a name="l00136"></a>00136 
<a name="l00137"></a>00137         <span class="keywordflow">while</span>(!q.empty()) {
<a name="l00138"></a>00138 
<a name="l00139"></a>00139                 <span class="comment">// Read the next entry on the queue</span>
<a name="l00140"></a>00140                 i=q.front();q.pop();
<a name="l00141"></a>00141                 j=q.front();q.pop();
<a name="l00142"></a>00142                 k=q.front();q.pop();
<a name="l00143"></a>00143 
<a name="l00144"></a>00144                 <span class="comment">// Check intersection of this image</span>
<a name="l00145"></a>00145                 <span class="keywordflow">if</span>(<a class="code" href="classvoro_1_1unitcell.html#acc445d2a964296dfca277e03fdc63ebd">intersects_image</a>(i,j,k,vol)) {
<a name="l00146"></a>00146 
<a name="l00147"></a>00147                         <span class="comment">// Add this entry to the output vectors</span>
<a name="l00148"></a>00148                         vi.push_back(i);
<a name="l00149"></a>00149                         vi.push_back(j);
<a name="l00150"></a>00150                         vi.push_back(k);
<a name="l00151"></a>00151                         vd.push_back(vol);
<a name="l00152"></a>00152 
<a name="l00153"></a>00153                         <span class="comment">// Add neighbors to the queue if they have not been</span>
<a name="l00154"></a>00154                         <span class="comment">// tested</span>
<a name="l00155"></a>00155                         ap=ac+i+ms2*(j+ms2*k);
<a name="l00156"></a>00156                         <span class="keywordflow">if</span>(k&gt;-max_unit_voro_shells&amp;&amp;*(ap-ms2*ms2)) {q.push(i);q.push(j);q.push(k-1);*(ap-ms2*ms2)=<span class="keyword">false</span>;}
<a name="l00157"></a>00157                         <span class="keywordflow">if</span>(j&gt;-max_unit_voro_shells&amp;&amp;*(ap-ms2)) {q.push(i);q.push(j-1);q.push(k);*(ap-ms2)=<span class="keyword">false</span>;}
<a name="l00158"></a>00158                         <span class="keywordflow">if</span>(i&gt;-max_unit_voro_shells&amp;&amp;*(ap-1)) {q.push(i-1);q.push(j);q.push(k);*(ap-1)=<span class="keyword">false</span>;}
<a name="l00159"></a>00159                         <span class="keywordflow">if</span>(i&lt;max_unit_voro_shells&amp;&amp;*(ap+1)) {q.push(i+1);q.push(j);q.push(k);*(ap+1)=<span class="keyword">false</span>;}
<a name="l00160"></a>00160                         <span class="keywordflow">if</span>(j&lt;max_unit_voro_shells&amp;&amp;*(ap+ms2)) {q.push(i);q.push(j+1);q.push(k);*(ap+ms2)=<span class="keyword">false</span>;}
<a name="l00161"></a>00161                         <span class="keywordflow">if</span>(k&lt;max_unit_voro_shells&amp;&amp;*(ap+ms2*ms2)) {q.push(i);q.push(j);q.push(k+1);*(ap+ms2*ms2)=<span class="keyword">false</span>;}
<a name="l00162"></a>00162                 }
<a name="l00163"></a>00163         }
<a name="l00164"></a>00164 
<a name="l00165"></a>00165         <span class="comment">// Remove mask memory</span>
<a name="l00166"></a>00166         <span class="keyword">delete</span> [] a;
<a name="l00167"></a>00167 }
<a name="l00168"></a>00168 <span class="comment"></span>
<a name="l00169"></a>00169 <span class="comment">/** Tests to see if a shell of periodic images could possibly cut the periodic</span>
<a name="l00170"></a>00170 <span class="comment"> * unit cell.</span>
<a name="l00171"></a>00171 <span class="comment"> * \param[in] l the index of the shell to consider.</span>
<a name="l00172"></a>00172 <span class="comment"> * \return True if a point in the shell cuts the cell, false otherwise. */</span>
<a name="l00173"></a>00173 <span class="keywordtype">bool</span> unitcell::unit_voro_intersect(<span class="keywordtype">int</span> l) {
<a name="l00174"></a>00174         <span class="keywordtype">int</span> i,j;
<a name="l00175"></a>00175         <span class="keywordflow">if</span>(unit_voro_test(l,0,0)) <span class="keywordflow">return</span> <span class="keyword">true</span>;
<a name="l00176"></a>00176         <span class="keywordflow">for</span>(i=1;i&lt;l;i++) {
<a name="l00177"></a>00177                 <span class="keywordflow">if</span>(unit_voro_test(l,i,0)) <span class="keywordflow">return</span> <span class="keyword">true</span>;
<a name="l00178"></a>00178                 <span class="keywordflow">if</span>(unit_voro_test(-l,i,0)) <span class="keywordflow">return</span> <span class="keyword">true</span>;
<a name="l00179"></a>00179         }
<a name="l00180"></a>00180         <span class="keywordflow">for</span>(i=-l;i&lt;=l;i++) <span class="keywordflow">if</span>(unit_voro_test(i,l,0)) <span class="keywordflow">return</span> <span class="keyword">true</span>;
<a name="l00181"></a>00181         <span class="keywordflow">for</span>(i=1;i&lt;l;i++) <span class="keywordflow">for</span>(j=-l+1;j&lt;=l;j++) {
<a name="l00182"></a>00182                 <span class="keywordflow">if</span>(unit_voro_test(l,j,i)) <span class="keywordflow">return</span> <span class="keyword">true</span>;
<a name="l00183"></a>00183                 <span class="keywordflow">if</span>(unit_voro_test(-j,l,i)) <span class="keywordflow">return</span> <span class="keyword">true</span>;
<a name="l00184"></a>00184                 <span class="keywordflow">if</span>(unit_voro_test(-l,-j,i)) <span class="keywordflow">return</span> <span class="keyword">true</span>;
<a name="l00185"></a>00185                 <span class="keywordflow">if</span>(unit_voro_test(j,-l,i)) <span class="keywordflow">return</span> <span class="keyword">true</span>;
<a name="l00186"></a>00186         }
<a name="l00187"></a>00187         <span class="keywordflow">for</span>(i=-l;i&lt;=l;i++) <span class="keywordflow">for</span>(j=-l;j&lt;=l;j++) <span class="keywordflow">if</span>(unit_voro_test(i,j,l)) <span class="keywordflow">return</span> <span class="keyword">true</span>;
<a name="l00188"></a>00188         <span class="keywordflow">return</span> <span class="keyword">false</span>;
<a name="l00189"></a>00189 }
<a name="l00190"></a>00190 <span class="comment"></span>
<a name="l00191"></a>00191 <span class="comment">/** Tests to see if a plane cut from a particular periodic image will cut th</span>
<a name="l00192"></a>00192 <span class="comment"> * unit Voronoi cell.</span>
<a name="l00193"></a>00193 <span class="comment"> * \param[in] (i,j,k) the index of the periodic image to consider.</span>
<a name="l00194"></a>00194 <span class="comment"> * \return True if the image cuts the cell, false otherwise. */</span>
<a name="l00195"></a>00195 <span class="keyword">inline</span> <span class="keywordtype">bool</span> unitcell::unit_voro_test(<span class="keywordtype">int</span> i,<span class="keywordtype">int</span> j,<span class="keywordtype">int</span> k) {
<a name="l00196"></a>00196         <span class="keywordtype">double</span> x=i*<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+j*<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>+k*<a class="code" href="classvoro_1_1unitcell.html#a000adf05e666f0fb4f57dbe8ca5b51db">bxz</a>,y=j*<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>+k*<a class="code" href="classvoro_1_1unitcell.html#a434c1a7cd8ba84d4ad09845a1a0ddb5b">byz</a>,z=k*<a class="code" href="classvoro_1_1unitcell.html#a30bf0a4a6bff808301aa424683eb0481">bz</a>;
<a name="l00197"></a>00197         <span class="keywordtype">double</span> rsq=x*x+y*y+z*z;
<a name="l00198"></a>00198         <span class="keywordflow">return</span> <a class="code" href="classvoro_1_1unitcell.html#a0844f70778cc7b19024ed28a2c7ad02c">unit_voro</a>.plane_intersects(x,y,z,rsq);
<a name="l00199"></a>00199 }
<a name="l00200"></a>00200 <span class="comment"></span>
<a name="l00201"></a>00201 <span class="comment">/** Draws the periodic domain in gnuplot format.</span>
<a name="l00202"></a>00202 <span class="comment"> * \param[in] fp the file handle to write to. */</span>
<a name="l00203"></a><a class="code" href="classvoro_1_1unitcell.html#a9379963d63388ba4075c620ddfdd326f">00203</a> <span class="keywordtype">void</span> <a class="code" href="classvoro_1_1unitcell.html#a630c635fd21cdd4e50bab118a4dadb13">unitcell::draw_domain_gnuplot</a>(FILE *fp) {
<a name="l00204"></a>00204         fprintf(fp,<span class="stringliteral">&quot;0 0 0\n%g 0 0\n%g %g 0\n%g %g 0\n&quot;</span>,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>,<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>);
<a name="l00205"></a>00205         fprintf(fp,<span class="stringliteral">&quot;%g %g %g\n%g %g %g\n%g %g %g\n%g %g %g\n&quot;</span>,<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>+bxz,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>+byz,bz,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>+bxz,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>+byz,bz,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+bxz,byz,bz,bxz,byz,bz);
<a name="l00206"></a>00206         fprintf(fp,<span class="stringliteral">&quot;0 0 0\n%g %g 0\n\n%g %g %g\n%g %g %g\n\n&quot;</span>,<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>,bxz,byz,bz,<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>+bxz,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>+byz,bz);
<a name="l00207"></a>00207         fprintf(fp,<span class="stringliteral">&quot;%g 0 0\n%g %g %g\n\n%g %g 0\n%g %g %g\n\n&quot;</span>,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+bxz,byz,bz,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>+bxz,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>+byz,bz);
<a name="l00208"></a>00208 }
<a name="l00209"></a>00209 <span class="comment"></span>
<a name="l00210"></a>00210 <span class="comment">/** Draws the periodic domain in POV-Ray format.</span>
<a name="l00211"></a>00211 <span class="comment"> * \param[in] fp the file handle to write to. */</span>
<a name="l00212"></a><a class="code" href="classvoro_1_1unitcell.html#a04c955e90a9478e17fb2a0c5df15a498">00212</a> <span class="keywordtype">void</span> <a class="code" href="classvoro_1_1unitcell.html#af476fca3dcc658abe2796aa9613f67bc">unitcell::draw_domain_pov</a>(FILE *fp) {
<a name="l00213"></a>00213         fprintf(fp,<span class="stringliteral">&quot;cylinder{0,0,0&gt;,&lt;%g,0,0&gt;,rr}\n&quot;</span>
<a name="l00214"></a>00214                    <span class="stringliteral">&quot;cylinder{&lt;%g,%g,0&gt;,&lt;%g,%g,0&gt;,rr}\n&quot;</span>,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>,<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>);
<a name="l00215"></a>00215         fprintf(fp,<span class="stringliteral">&quot;cylinder{&lt;%g,%g,%g&gt;,&lt;%g,%g,%g&gt;,rr}\n&quot;</span>
<a name="l00216"></a>00216                    <span class="stringliteral">&quot;cylinder{&lt;%g,%g,%g&gt;,&lt;%g,%g,%g&gt;,rr}\n&quot;</span>,bxz,byz,bz,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+bxz,byz,bz,<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>+bxz,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>+byz,bz,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>+bxz,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>+byz,bz);
<a name="l00217"></a>00217         fprintf(fp,<span class="stringliteral">&quot;cylinder{&lt;0,0,0&gt;,&lt;%g,%g,0&gt;,rr}\n&quot;</span>
<a name="l00218"></a>00218                    <span class="stringliteral">&quot;cylinder{&lt;%g,0,0&gt;,&lt;%g,%g,0&gt;,rr}\n&quot;</span>,<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>);
<a name="l00219"></a>00219         fprintf(fp,<span class="stringliteral">&quot;cylinder{&lt;%g,%g,%g&gt;,&lt;%g,%g,%g&gt;,rr}\n&quot;</span>
<a name="l00220"></a>00220                    <span class="stringliteral">&quot;cylinder{&lt;%g,%g,%g&gt;,&lt;%g,%g,%g&gt;,rr}\n&quot;</span>,bxz,byz,bz,<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>+bxz,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>+byz,bz,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+bxz,byz,bz,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>+bxz,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>+byz,bz);
<a name="l00221"></a>00221         fprintf(fp,<span class="stringliteral">&quot;cylinder{&lt;0,0,0&gt;,&lt;%g,%g,%g&gt;,rr}\n&quot;</span>
<a name="l00222"></a>00222                    <span class="stringliteral">&quot;cylinder{&lt;%g,0,0&gt;,&lt;%g,%g,%g&gt;,rr}\n&quot;</span>,bxz,byz,bz,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+bxz,byz,bz);
<a name="l00223"></a>00223         fprintf(fp,<span class="stringliteral">&quot;cylinder{&lt;%g,%g,0&gt;,&lt;%g,%g,%g&gt;,rr}\n&quot;</span>
<a name="l00224"></a>00224                    <span class="stringliteral">&quot;cylinder{&lt;%g,%g,0&gt;,&lt;%g,%g,%g&gt;,rr}\n&quot;</span>,<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>,<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>+bxz,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>+byz,bz,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>+bxz,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>+byz,bz);
<a name="l00225"></a>00225         fprintf(fp,<span class="stringliteral">&quot;sphere{&lt;0,0,0&gt;,rr}\nsphere{&lt;%g,0,0&gt;,rr}\n&quot;</span>
<a name="l00226"></a>00226                    <span class="stringliteral">&quot;sphere{&lt;%g,%g,0&gt;,rr}\nsphere{&lt;%g,%g,0&gt;,rr}\n&quot;</span>,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>,<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>);
<a name="l00227"></a>00227         fprintf(fp,<span class="stringliteral">&quot;sphere{&lt;%g,%g,%g&gt;,rr}\nsphere{&lt;%g,%g,%g&gt;,rr}\n&quot;</span>
<a name="l00228"></a>00228                    <span class="stringliteral">&quot;sphere{&lt;%g,%g,%g&gt;,rr}\nsphere{&lt;%g,%g,%g&gt;,rr}\n&quot;</span>,bxz,byz,bz,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+bxz,byz,bz,<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>+bxz,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>+byz,bz,<a class="code" href="classvoro_1_1unitcell.html#aba3ac9d30b539fb0017799555199c2df">bx</a>+<a class="code" href="classvoro_1_1unitcell.html#a572a7b91539af3c0707b6caccbf15f35">bxy</a>+bxz,<a class="code" href="classvoro_1_1unitcell.html#abaca9a0c853b9a53f164a2fbee783f83">by</a>+byz,bz);
<a name="l00229"></a>00229 }
<a name="l00230"></a>00230 
<a name="l00231"></a>00231 }
</pre></div></div>
</div>


<hr class="footer"/><address class="footer"><small>
Generated on Fri Sep 23 2011 22:49:06 for Voro++ by &#160;<a href="http://www.doxygen.org/index.html">
<img class="footer" src="doxygen.png" alt="doxygen"/>
</a> 1.7.5.1
</small></address>

</body>
</html>
